function stats = BondRetStats_nominal(parms,gesol,prob_NxYz,Phi_NxYz_1y,grid_Y,stats_in,FLAG_RUN_PARALLEL)

    if nargin<7
        FLAG_RUN_PARALLEL = false;
    end

    stats = stats_in;
    
    mats = 2:5;

    [NN, Nx, NY, Nz] = size(prob_NxYz);

    idx_1y = find(gesol.nom_bonds.maturities==12);
    p1_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_1y));
    Y_NxYz = permute(repmat(grid_Y(:),[1 NN Nx Nz]),[2 3 1 4]);
    p_1y_NxYz = permute(repmat(p1_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        - gesol.nom_bonds.P_C_a(idx_1y)*Y_NxYz;

    bond_risk_prem_uncond = zeros(size(mats));
    bond_exret_uncond_vol = zeros(size(mats));
    bond_ret_uncond_vol = zeros(size(mats));
    
    if FLAG_RUN_PARALLEL
        
        % prepare inputs
        term0a_list = cell(numel(mats),1);
        term0b_list = cell(numel(mats),1);
        termt_list = cell(numel(mats),1);
        inflation_process = gesol.inflation_process;
        
        for n = 1:numel(mats)
            
            mat = mats(n);
            idx_buy = find(gesol.nom_bonds.maturities==12*mat);
            idx_sell = find(gesol.nom_bonds.maturities==12*(mat-1));

            pbuy_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_buy));
            psell_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));
            
            % p_sell - p_buy + p_1y
            term0.fun_Nxz = -pbuy_Nxz + p1_Nxz;
            term0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_1y);
            term0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_b(idx_1y);
            term0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);
            term0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            termt.fun_Nxz = psell_Nxz;
            termt.const = -gesol.nom_bonds.P_B_a(idx_sell);
            termt.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            termt.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);
            termt.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            
            term0a_list{n} = term0;
            termt_list{n} = termt;
            
            % p_sell - p_buy
            term0.fun_Nxz = -pbuy_Nxz;
            term0.const = gesol.nom_bonds.P_B_a(idx_buy);
            term0.slope_X = gesol.nom_bonds.P_B_b(idx_buy);
            term0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy);
            term0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy);
            
            term0b_list{n} = term0;
            
        end
        
        parfor n = 1:numel(mats)

            % p_sell - p_buy + p_1y
            [M,VAR] = generic_nominal_moments(inflation_process,...
                term0a_list{n},12,termt_list{n},prob_NxYz, Phi_NxYz_1y, grid_Y);
            bond_risk_prem_uncond(n) = M;
            bond_exret_uncond_vol(n) = sqrt(VAR);

            % p_sell - p_buy
            [M,VAR] = generic_nominal_moments(inflation_process,...
                term0b_list{n},12,termt_list{n},prob_NxYz, Phi_NxYz_1y, grid_Y);
            bond_ret_uncond_vol(n) = sqrt(VAR);

        end
        
    else

        for n = 1:numel(mats)

            mat = mats(n);
            idx_buy = find(gesol.nom_bonds.maturities==12*mat);
            idx_sell = find(gesol.nom_bonds.maturities==12*(mat-1));

            pbuy_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_buy));
            psell_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));

            % p_sell - p_buy + p_1y
            term0.fun_Nxz = -pbuy_Nxz + p1_Nxz;
            term0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_1y);
            term0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_b(idx_1y);
            term0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);
            term0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            termt.fun_Nxz = psell_Nxz;
            termt.const = -gesol.nom_bonds.P_B_a(idx_sell);
            termt.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            termt.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);
            termt.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            [M,VAR] = generic_nominal_moments(gesol.inflation_process,...
                term0,12,termt,prob_NxYz, Phi_NxYz_1y, grid_Y);
            bond_risk_prem_uncond(n) = M;
            bond_exret_uncond_vol(n) = sqrt(VAR);

            % p_sell - p_buy
            term0.fun_Nxz = -pbuy_Nxz;
            term0.const = gesol.nom_bonds.P_B_a(idx_buy);
            term0.slope_X = gesol.nom_bonds.P_B_b(idx_buy);
            term0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy);
            term0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy);
            [M,VAR] = generic_nominal_moments(gesol.inflation_process,...
                term0,12,termt,prob_NxYz, Phi_NxYz_1y, grid_Y);
            bond_ret_uncond_vol(n) = sqrt(VAR);

        end
    
    end

    stats.bond_risk_prem_nom_uncond = bond_risk_prem_uncond;
    
    % conditional bond risk premia, conditional on z
    bond_risk_prem_cond = zeros(numel(mats),parms.Nz);
    
    prob_NxY_given_z = zeros(parms.NN,parms.Nx,NY,parms.Nz);
    for iz = 1:parms.Nz
        prob_NxY_given_z(:,:,:,iz) = prob_NxYz(:,:,:,iz)/sum(prob_NxYz(:,:,:,iz),'all');
    end

    mean_X = gesol.inflation_process.mu_pi/(1 - gesol.inflation_process.rho_pi);
    mean_p_1y_Xe = -gesol.nom_bonds.P_B_a(idx_1y) - gesol.nom_bonds.P_B_b(idx_1y)*mean_X;
    
    for n = 1:numel(mats)
        mat = mats(n);
        idx_buy = find(gesol.nom_bonds.maturities==mat*12);
        idx_sell = find(gesol.nom_bonds.maturities==(mat-1)*12);
        p_buy_NxYz = permute(repmat(log(gesol.nom_bonds.P_A(:,:,:,idx_buy)),[1 1 1 NY]),[1 2 4 3]) ...
            - gesol.nom_bonds.P_C_a(idx_buy)*Y_NxYz;
        mean_p_buy_Xe = -gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_b(idx_buy)*mean_X;
        p_sell_NxYz = permute(repmat(log(gesol.nom_bonds.P_A(:,:,:,idx_sell)),[1 1 1 NY]),[1 2 4 3]) ...
            - gesol.nom_bonds.P_C_a(idx_sell)*Y_NxYz;
        mean_p_sell_Xe = -gesol.nom_bonds.P_B_a(idx_sell) - gesol.nom_bonds.P_B_b(idx_sell)*mean_X;
        E_p_sell_1y_NxYz = reshape(Phi_NxYz_1y*p_sell_NxYz(:),[parms.NN parms.Nx NY parms.Nz]);
        if FLAG_RUN_PARALLEL
            parfor iz = 1:parms.Nz
                bond_risk_prem_cond(n,iz) = sum(prob_NxY_given_z(:,:,:,iz).*(E_p_sell_1y_NxYz(:,:,:,iz) ...
                    - p_buy_NxYz(:,:,:,iz) + p_1y_NxYz(:,:,:,iz)),'all') ...
                    + mean_p_sell_Xe - mean_p_buy_Xe + mean_p_1y_Xe;
            end
        else
            for iz = 1:parms.Nz
                bond_risk_prem_cond(n,iz) = sum(prob_NxY_given_z(:,:,:,iz).*(E_p_sell_1y_NxYz(:,:,:,iz) ...
                    - p_buy_NxYz(:,:,:,iz) + p_1y_NxYz(:,:,:,iz)),'all') ...
                    + mean_p_sell_Xe - mean_p_buy_Xe + mean_p_1y_Xe;
            end
        end
    end
    
    for iz = 1:parms.Nz
        varname = ['bond_risk_prem_nom_z',num2str(iz)];
        stats.(varname) = reshape(bond_risk_prem_cond(:,iz),1,[]);
    end
    
    stats.bond_ret_nom_uncond_vol = bond_ret_uncond_vol;
    stats.bond_exret_nom_uncond_vol = bond_exret_uncond_vol;
    
end

%% old save
% 
%     inflation_process = gesol.inflation_process;
% 
%     mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
%     var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi ...
%         + inflation_process.theta_pi^2)*(inflation_process.sigma_pi^2)...
%         /(1 - inflation_process.rho_pi^2);
%     var_eps = inflation_process.sigma_pi^2;
%     cov_Xe = var_eps;
% 
%     cov_X0_X12 = (inflation_process.rho_pi^12)*var_X ...
%         + (inflation_process.rho_pi^11)*inflation_process.theta_pi*cov_Xe;
%     cov_eps0_X12 = (inflation_process.rho_pi^12)*cov_Xe ...
%         + (inflation_process.rho_pi^11)*inflation_process.theta_pi*var_eps;
% 
%     Y_NxYz = permute(repmat(grid_Y(:),[1 parms.NN parms.Nx parms.Nz]),[2 3 1 4]);
%     NY = numel(grid_Y);
% 
%     % unconditional bond risk premia
%     mats = 2:5; % years
%     bond_risk_prem_uncond = zeros(size(mats));
% 
%     % unconditional volatilities
%     bond_ret_uncond_vol = zeros(size(mats));
%     bond_exret_uncond_vol = zeros(size(mats));
% 
%     idx_1y = find(gesol.nom_bonds.maturities==12);
%     p_1y_NxYz = permute(repmat(log(gesol.nom_bonds.prices(:,:,:,idx_1y)),[1 1 1 NY]),[1 2 4 3]) ...
%         - gesol.nom_bonds.P_C_a(idx_1y)*Y_NxYz;
%     mean_p_1y_NxYz = sum(prob_NxYz.*p_1y_NxYz,'all');
%     var_p_1y_NxYz = sum(prob_NxYz.*(p_1y_NxYz.^2),'all') - mean_p_1y_NxYz^2;
%     mean_p_1y_Xe = -gesol.nom_bonds.P_B_a(idx_1y) - gesol.nom_bonds.P_B_b(idx_1y)*mean_X;
%     var_p_1y_Xe = (gesol.nom_bonds.P_B_b(idx_1y)^2)*var_X + (gesol.nom_bonds.P_B_c(idx_1y)^2)*var_eps ...
%         + 2*gesol.nom_bonds.P_B_b(idx_1y)*gesol.nom_bonds.P_B_c(idx_1y)*cov_Xe;
% 
%     for n = 1:numel(mats)
%         % find indices
%         mat = mats(n);
%         idx_buy = find(gesol.nom_bonds.maturities==mat*12);
%         idx_sell = find(gesol.nom_bonds.maturities==(mat-1)*12);
%         % p_buy and p_sell
%         p_buy_NxYz = permute(repmat(log(gesol.nom_bonds.prices(:,:,:,idx_buy)),[1 1 1 NY]),[1 2 4 3]) ...
%             - gesol.nom_bonds.P_C_a(idx_buy)*Y_NxYz;
%         p_sell_NxYz = permute(repmat(log(gesol.nom_bonds.prices(:,:,:,idx_sell)),[1 1 1 NY]),[1 2 4 3]) ...
%             - gesol.nom_bonds.P_C_a(idx_sell)*Y_NxYz;
%         % E[p] and Var[p], NxYz component
%         mean_p_buy_NxYz = sum(prob_NxYz.*p_buy_NxYz,'all');
%         var_p_buy_NxYz = sum(prob_NxYz.*(p_buy_NxYz.^2),'all') - mean_p_buy_NxYz^2;
%         mean_p_sell_NxYz = sum(prob_NxYz.*p_sell_NxYz,'all');
%         var_p_sell_NxYz = sum(prob_NxYz.*(p_sell_NxYz.^2),'all') - mean_p_sell_NxYz^2;
%         % E[p] and Var[p], Xe component
%         mean_p_buy_Xe = -gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_b(idx_buy)*mean_X;
%         var_p_buy_Xe = (gesol.nom_bonds.P_B_b(idx_buy)^2)*var_X + (gesol.nom_bonds.P_B_c(idx_buy)^2)*var_eps ...
%             + 2*gesol.nom_bonds.P_B_b(idx_buy)*gesol.nom_bonds.P_B_c(idx_buy)*cov_Xe;
%         mean_p_sell_Xe = -gesol.nom_bonds.P_B_a(idx_sell) - gesol.nom_bonds.P_B_b(idx_sell)*mean_X;
%         var_p_sell_Xe = (gesol.nom_bonds.P_B_b(idx_sell)^2)*var_X + (gesol.nom_bonds.P_B_c(idx_sell)^2)*var_eps ...
%             + 2*gesol.nom_bonds.P_B_b(idx_sell)*gesol.nom_bonds.P_B_c(idx_sell)*cov_Xe;
%         % Var[p_sell-p_buy]
%         var_ret_NxYz = var_p_sell_NxYz + var_p_buy_NxYz ...
%             - 2*(sum(prob_NxYz(:).*p_buy_NxYz(:).*(Phi_NxYz_1y*p_sell_NxYz(:))) - mean_p_sell_NxYz*mean_p_buy_NxYz);
%         cov_p_buy_p_sell_Xe = gesol.nom_bonds.P_B_b(idx_sell)*gesol.nom_bonds.P_B_b(idx_buy)*cov_X0_X12 ...
%             + gesol.nom_bonds.P_B_b(idx_sell)*gesol.nom_bonds.P_B_c(idx_buy)*cov_eps0_X12;
%         var_ret_Xe = var_p_sell_Xe + var_p_buy_Xe - 2*cov_p_buy_p_sell_Xe;
%         % Var[p_sell-p_buy+p_1y]
%         mean_p_buy_p_1y_NxYz = -mean_p_buy_NxYz + mean_p_1y_NxYz;
%         var_p_buy_p_1y_NxYz = sum(prob_NxYz.*((-p_buy_NxYz + p_1y_NxYz).^2),'all') - mean_p_buy_p_1y_NxYz^2;
%         var_exret_NxYz = var_p_sell_NxYz + var_p_buy_p_1y_NxYz ...
%             + 2*(sum(prob_NxYz(:).*(-p_buy_NxYz(:) + p_1y_NxYz(:)).*(Phi_NxYz_1y*p_sell_NxYz(:))) - mean_p_buy_p_1y_NxYz*mean_p_sell_NxYz);
%         var_p_buy_p_1y_Xe = ((gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_b(idx_1y))^2)*var_X ...
%             + ((gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y))^2)*var_eps ...
%             + 2*(gesol.nom_bonds.P_B_b(idx_buy)-gesol.nom_bonds.P_B_b(idx_1y))...
%             *(gesol.nom_bonds.P_B_c(idx_buy)-gesol.nom_bonds.P_B_c(idx_1y))*cov_Xe;
%         var_exret_Xe = var_p_sell_Xe + var_p_buy_p_1y_Xe + ;
%         % bond risk premium
%         bond_risk_prem_uncond(n) = mean_p_sell_NxYz - mean_p_buy_NxYz + mean_p_1y_NxYz ...
%             + mean_p_sell_Xe - mean_p_buy_Xe + mean_p_1y_Xe;
%         % return volatility
%         bond_ret_uncond_vol(n) = var_ret_NxYz + var_ret_Xe;
%         % excess return volatility
%         bond_exret_uncond_vol(n) = var_exret_NxYz + var_exret_Xe;
%     end
%     stats.bond_risk_prem_nom_uncond = bond_risk_prem_uncond;